*****************************************************************
* Replication directory for                                   ***
* Prime locations                                             ***
* by Gabriel M. Ahlfeldt, Thilo N.H. Albers, Kristian Behrens ***
* Published in American Economic Review: Insights             ***
*****************************************************************
* 01/2025
* Stata
version 17.0

* This do file generates internal and external distance from prime location boundaries in Global Cities
* These distances will be used in the BDD plots that show how emplyoment densities change across the prime location boundary 

* Generate temporary working directory
	qui capture mkdir "$temp/TEMPDISTBORDER" // Generate temporary working directory

* Read metro identifier for loop 	
	use "$temp/GIS_using/PL125_grids", clear
	levelsof metro_id	, local(MID)

* Start loop over metros	
	foreach m_ID of local MID {
	use "$temp/GIS_using/PL125_grids", clear
	keep if metro_id==`m_ID'
	* identify neighbor cells with Euclidian distance: a border cell is a cell with less than at neighbors 
		preserve
		keep square_id grid_x grid_y 
		rename grid_x x
		rename grid_y y
		tempfile points
		save `points'
		rename * *0
		cross using `points'
		gen double d=sqrt((x-x0)^2+(y-y0)^2)
		*everything below 354 meters is an adjacent cell
		drop if d==0
		gen neighbor_i=1 if d<354
		egen neighbors=sum(neighbor_i), by(square_id)
		gen border_cell=1 if neighbors!=8
		gen border_cell_corner=1 if neighbors==7
		replace border_cell=0 if border_cell==.
		keep square_id border_cell border_cell_corner
		duplicates drop square_id, force 
		tempfile neighbors_ids
		save `neighbors_ids'
		restore
		merge 1:1 square_id using `neighbors_ids', nogen
		* Calculate internal distances 
		preserve 
		keep if border_cell==1
		keep square_id lat lon 
		rename square_id nborid
		rename lat nborlat
		rename lon nborlon
		tempfile nbordata_border
		save  `nbordata_border' 
		restore 
		geonear  square_id lat lon using `nbordata_border' , genstub(borderdist_) n(nborid nborlat nborlon) near(1) 
		* Calculate correction factor for converting planar to geodesic distance
		gen correction_factor=cos(2*atan(exp([grid_y^2^0.5]/6371000))-_pi/2)
		* add the geodeisc equivalent of 125m planar distance to all cells to account for difference between centroid and distance to actual border of cell  
		gen border_dist_internal=-(km_to_borderdist_+correction_factor*.125)
		*replace for those that are only corner 
		replace border_dist_internal=-(correction_factor*(sqrt(.125^2+.125^2))) if  border_cell_corner==1
		* clear up variables 
		drop borderdist_ km_to_borderdist_         correction_factor
		*** calculate external distances
		* join full city grid
		merge 1:1 grid_x grid_y  using    "$temp/125CitiesPLs/grid125_PL_output.dta"
		drop if metro_id!=`m_ID'
		* calulate correction factor for geodesic vs. planar
		gen correction_factor=cos(2*atan(exp([grid_y^2^0.5]/6371000))-_pi/2)
		drop _merge 
		geonear  square_id lat lon using `nbordata_border'  , genstub(BD_external) n(nborid nborlat nborlon) near(1)  
		sort km_to_BD_external
		*substract   125 meters as we are calculating centroid to centroid distances here
		gen border_dist_external=km_to_BD_external-(correction_factor*.125) if PL==0 
		gen distance_toPLborder=border_dist_internal if PL==1
		replace distance_toPLborder=border_dist_external if PL==0
		keep grid_x grid_y square_id distance_toPLborder
 		save "$temp/TEMPDISTBORDER/DIST_`m_ID'", replace 					  
	}
 * Loop ends
 
 * Generate whole grid data and export
	clear
	foreach m_ID of local MID {	// Append all metro data sets
	append using "$temp/TEMPDISTBORDER/DIST_`m_ID'"
	}
	merge 1:1 grid_x grid_y square_id using  "$temp/125CitiesPLs/grid125_PL_output.dta", nogen  // Merge grid data
 
 * Clean up and save
	preserve
		filelist, dir("$temp/TEMPDISTBORDER")
		levelsof filename, local(files)
		foreach f of local files{
			erase "$temp/TEMPDISTBORDER/`f'"
		}
	restore 
 	label var distance_toPLborder "distance to border of PL (internal distances negative)"
	save  "$temp/grid125_PL_distances", replace

* Script ends	